A model for the compaction of granular media 
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We introduce a lattice model, in which frustration plays 
a crucial role, to describe relaxation properties of granular 
media. We show Monte Carlo results for compaction in the 
presence of vibrations and gravity, which compare well with 
experimental data. 

Despite their importance for industrial applications re- 
laxation phenomena in non-thermal disordered systems 
as granular media, have just recently begun to be stud- 
ied systematically. A common and simple experiment in 
this context, is the compaction of sand. When a box 
filled with loose packed sand is shaken at low amplitude, 
density visibly increases. If in addition the density goes 
beyond a definite threshold the mechanical properties of 
sand abruptly change and the granular structure can- 
not be sheared any longer without a volume increase. 
This phenomenon, very important in practical applica- 
tions Q, was observed by Reynolds Q and is referred 
to as the "Reynolds" or "dilatancy" transition. 

For a given macroscopic parameter as density a gran- 
ular packing can be in a huge number of different mi- 
croscopical states. In order to describe this situation 
concepts from statistical mechanics have been introduced 
Relations to spin glasses (SG) have been suggested 
several years ago (see references in ||). In fact a char- 
acteristic of SG is their non trivial phase space which 
gives rise to its complex static and dynamic behavior. 
The phase space structure of SG is due to the presence 
of quenched disorder and frustration. Strictly speaking 
quenched disorder is not present in granular media but 
there are effects of "geometric frustration", known also 
from hard sphere systems. This kind of frustration is 
generated by the steric constraints imposed by the hard 
core repulsion of neighboring grains and the subsequent 
interlocking which leads to non local cooperative macro- 
scopic rearrangements. Recently the analogies between 
an intrinsically frustrated system like frustrated percola- 
tion B and phase transitions in a granular packing have 
been outlined JtJ. 

In this paper we present computer simulations of a sim- 
ple frustrated Ising lattice gas model, subject to gravity 
following a diffusion like Monte Carlo dynamics. The 
particles in this system are characterized by internal de- 
grees of freedom which describe their orientation or other 
local sterical properties of the grains. This model with- 
out gravity shows complex behavior similar to the one 
observed in glass forming liquids and spin glasses M . We 



will show how the density of our lattice gas is strongly 
dependent on the duration and the amplitude of simply 
implemented vibrations. Our data reproduce the loga- 
rithmic relaxation behavior found in real experiments in 
a sequence of taps and offer the possibility to make new 
predictions also for single tap processes. Our data also 
reproduce the distribution of forces at the bottom of the 
system as found in real experiments. A relation appears 
between the SG transition, signaled by the vanishing of 
macroscopic self-diffusion, and the Reynolds transition in 
granular systems. 

We consider a system of particles which move on a 
square lattice whose bonds are characterized by quenched 
random numbers e.y = ±1. On site i we set n; = 1 if a 
particle is present and otherwise. The particles have 
an internal degree of freedom Si — ±1 and are subjected 
to the constraint that whenever two (i and j) are neigh- 
boring, their "spin" must satisfy the relation 
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i.e. they have to fit the local "geometrical" structure. 
When the density of particles is high enough they can 
feel the frustration that has been imposed by the choice 
of the 6ij . As a consequence, in resemblance to frustrated 
percolation Q , particles can never close a frustrated loop 
in the lattice leaving empty sites (see below). 

The physical origin of the bond variables Ey, is the 
geometrical frustration originated in granular systems by 
the actual shapes and arrangements of particles and the 
internal variables Si mimic local shapes and positions. 

We have studied this system when subject to "gravity" 
and "external vibrations" . The dynamics of our model 
consists in a random diffusion of particles on a square lat- 
tice tilted by 45° (see Fig.p in such a way as to preserve 
the above constraint. The particles attempt a move up- 
ward with probability Pi and downward with P\ (with 
Pi + P2 = 1). The move is made only if the internal 
degrees of freedom satisfy eq. (|f|). Similarly a spin flips 
with probability one if there is no violation of eq. (Q), 
and does not flip otherwise. In absence of vibrations, the 
effect of gravity imposes P2 — 0. When vibrations are 
switched on Pi becomes finite. The crucial parameter 
which controls the dynamics and the final density is the 
ratio x(t) = P2{t)/P\(t) which describes the amplitude 
of the vibration. 

It is possible to associate to this model a standard 
Hamiltonian formalism and establish a magnetic analogy, 
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based on the following definition: 

- H = y^^J(ejjSjSj - l)riinj + (2) 

(y> « 

where Sj = ±1 are spin variables, rij = 0, 1 occupancy 
variables and e,j = ±1 quenched interactions associated 
to the bonds of the lattice. It has been shown in mean 
field approximation |j) and numerically for finite dimen- 
sional systems ||, that Hamiltonian (||) exhibits a spin 
glass transition at high density (or low temperature). 

This Hamiltonian reduces in the p — > oo limit (all sites 
occupied) to the usual ±J Jsmj spin glass |To| . In the 
limit J — > oo, it describes a lattice gas in which frus- 
trated loops entirely filled with sites are forbidden be- 
cause along closed paths energetic reasons impose the 

quantity ioop( e v^3 ~ 1 ) to be zera In tnis umit 

a version of site frustrated percolation is recovered . 
When the particle number is fixed the configuration space 
of the system obtained in this last limit is the same as 
that of the frustrated lattice gas introduced at the begin- 
ning of this paragraph. 

We have studied the model introduced above in a 2D 
box with periodic boundary conditions along the x-axis 
and rigid walls at its bottom and top. 

After fixing the random quenched on the bonds, a 
random initial particle configuration is prepared by ran- 
domly inserting particles of given spin into the box from 
its top and then letting them fall down, with the de- 
scribed dynamics (P2 = 0), until the box is filled. To 
obtain an initial low density configuration we do not al- 
low particle spins to flip in this preparation process. The 
state prepared in this way has a density of about 0.520 
and corresponds to a random loose packing. 

We know experimentally that sand which is randomly 
poured into a box reaches higher density after shaking. 
In some experiments, the shaking process occurs in a se- 
quence of "taps". A tap is defined by its duration and 
its amplitude. After a sequence of taps the density de- 
cays logarithmically to a static limit (see 0| ) . We have 
studied the phenomena of density relaxation during a 
sequence of taps. In our MC simulation, each tap is a 
process in which vibrations are step like: x(t) = xq if 
t £ [0, t] while for t > r the system evolves subject only 
to gravity (Pi = 1 and P2 = 0, i.e. x(t) = 0) until 
it reaches a final "static" configuration JT^]. After each 
tap we have measured the bulk density of the system 
p(r, x ;t n ) defined as the mean density in the lower 25% 
of the box (t n is the n-th tap number) . Our results for 
density relaxation, in a box of size 30 x 60 averaged over 
32 different {e^ } configurations, are shown in the insert 
of Fig. |[ The behavior of p(r, x \t n ) is well fitted by 
the following logarithmic function in agreement with the 
experimental data (see |ll"| , |l3|| ): 

p(r, x a ;t) = p s - A Poo /[l + B ln(t/r + 1)] (3) 

In Fig. H we have collapsed our results for four different 
amplitudes as well as the experimental data for three 



different amplitudes on a single curve using eq. (0) and 
see that the agreement is very satisfactory. 

We have also simulated a single tapping process. In 
this case we have found that the relaxation, instead of 
being logarithmic as in the sequence of taps is well de- 
scribed by a "stretched exponential" . In the simula- 
tion we start at t = from a random loose packing 
configuration described before, then we introduce "vi- 
brations" in the interval t g [0, r] linearly decreasing 
the ratio x(t), x(t) = xq(1 — t/r) (with xq — 1). For 
t > t we put x(t) =0 and let the system evolve until it 
reaches a final "static" configuration. The final "static" 
bulk density p{j) monotonically increases with the vibra- 
tion time r asymptotically reaching a maximal density 
value p* ~ 0.78 when r — > 00 Q. 

During the dynamical process described above, we 
have recorded the time dependence of the mean bulk den- 
sity p(t,r). We find that the static limit is reached with 
a stretched "relaxation form" juj : 

p(t, t) = Ps {t) - A exp(-((t - t )/Tf) (4) 

Typical values of the parameters of eq. (Q) in our range 
of t are A 6 [0.15, 0.25], t € [-10 1 , 10 3 ], T e [10 2 , 10 4 ], 
j3 6 [2,4]. Note that the stretched exponential behavior 
only sets in after a time to, which can be very long if r 
is long. The relaxation processes found here are rather 
different from the logarithmic relaxation found in the se- 
quence of taps and could be investigated experimentally. 

The effect of compaction is clearly shown by the final 
density profile as a function of depth h, p(h, r), depicted 
in Fig^. In this case the box has a size 100 x 200 and the 
final states have been averaged over 32 to 512 different 
{ey} configurations (according the value of r). As sug- 
gested in ref. |l5j the density profile of granular media 
can be fitted using a generalized Fermi-Dirac distribu- 
tion. As shown in Fig.|3| the data from our model are well 
fitted by such a function for different values of r: 

p(h,T) = p s (T)[l-l/[l + exp((h-h (T))/s(T))}} . (5) 

To characterize a particle packing and its capability of 
internal rearrangement, we studied their self-diffusivity 
at fixed global particle density by setting x = 1. Specifi- 
cally we have studied the time dependence of the particle 
mean square displacement R 2 {t) — (j* J2i( r i(t)~ r i(®)) 2 ) 

A very interesting phenomenon is observed for densi- 
ties close to the maximal value p*: R 2 (t) shows devi- 
ations from the linear time dependence typical of stan- 
dard Brownian diffusive motion and presents an inflection 
point H . This signals the existence of two characteris- 
tic time regimes for particle motion (as already argued 
in p6|). From the long time behavior of R 2 (t) ~ Dt we 
extract the diffusion coefficient D(p), which goes to zero 
at about p*, signaling a localization transition in which 
particles are confined in local cages and the macroscopic 
diffusion-like processes are suppressed. This phenomenon 
may also be described in a different way: p* is the density 
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above which it becomes impossible to obtain a macro- 
scopic rearrangement of the particle positions without 
increasing the system volume, i.e. the density at which 
macroscopic shear in the system is impossible without 
dilatancy. This then seems to correspond to the quoted 
Reynolds transition in real granular media. 

The density p* coincides with the density at which the 
spin glass (SG) transition of Hamiltonian (||) (for J — > 
oo ) is located. This implies that at p* the SG correlation 
length £sg diverges, signaling the presence of collective 
behavior in the system. In SG this length can be detected 
by measuring the non linear susceptibility. In granular 
material the spin variables represent internal degrees of 
freedom and cannot be easily detected. 

The coincidence of the SG transition and the sup- 
pression of self-diffusivity suggest the equivalence of the 
Reynolds transition in granular media, the SG transi- 
tion in magnetic systems and the "ideal" glass transition 
in glass forming liquids [0J|]. Like in glass forming liq- 
uids the Reynolds transition does not show divergences 
in thermodynamic functions. 

The model introduced here is suited to study also other 
aspects of granular media. If a force is applied at the 
top of a granular system in a box, the local distribu- 
tion of forces v at the bottom follows an exponential law 
P(v) = a-exp(—c v). As suggested in |17H19|, it is possi- 
ble to introduce simplified models to describe the physics 
of forces in granular systems. In particular, a model has 
been proposed in JL8[ in which forces are represented by 
scalars. In our lattice model we apply the same approx- 
imation to study the force distributions in static config- 
urations of the system. We suppose that each present 
site (rii = 1) carries its own weight (equal to unity) and 
transmits the force Wi acting on it to its first left and 
right neighboring sites in the lower row. If its right (left) 
neighbor has a distance l r (If) from site i, the force con- 
tribution it receives from this site is equal to Wi-lij (l r +k) 
(wi ■ l r /(l r + If) respectively for the left site). We have 
calculated the force distribution P(v) at the bottom of 
our system as shown in Fig. ^. In agreement with the ex- 
perimental data and the result of the model introduced 
in , our data are well fitted by: 



P(v) = a ■ v b exp(—c v) 



(6) 



As noted in Ref. jig], the power law in front of the 
exponential would be very difficult to be detected exper- 
imentally since it affects the distribution for small value 
of v. 

In conclusion in this paper a frustrated Ising lattice gas 
has been introduced to describe different aspects of the 
phenomenology of granular systems, such as compaction 
in the presence of vibrations, logarithmic density relax- 
ation, and exponential force distribution. The results are 
in agreement with real experiments. The model is able 
to predict new results which arc amenable to experimen- 
tal observation, some of which have been reported here, 
while others are under investigation pM- The model 



which contains geometrical frustration as essential ingre- 
dient also shares features of spin glasses and glass forming 
liquids. 

Although we have reported here numerical results in 
2D we expect the same features also in 3D. 

We thank H. Jaeger and J. Knight for sending us their 
experimental data and IDRIS for computer time on Cray- 
T3D. 



[1] 

[2] 

[3] 



[7] 
[8] 
[9] 

[10] 

[11] 
[12] 



[13] 
[14] 
[15] 
[16] 
[17] 
[18] 

[19] 



Bashir Y.M. and Goddard J.D., J. Rheol. 35 (1991) 849; 
Vermeer P.A. and de Borst, Heron 29 (1984) 1. 
Reynolds O., Phil. Mag. S. 20 (1885) 469. 
H.M. Jaeger and S.R. Nagel, Science 255, 1523 (1992); 
H.M. Jaeger, S.R. Nagel and R.P. Behringer, Phys. To- 
day April 1996. 

Edwards S. F., J. Stat. Phys. 62 (1991) 889; Mehta A. 
and Edwards S. F., Physica A 157 (1989) 1091. 
H.J. Herrmann, J. Physique II 3, (1993) 427. 
Coniglio A., II Nuovo Cimento 16D (1994) 1027; Coniglio 
A., in Correlation and Connectivity, eds. H.E. Stanley 
and N. Ostrowsky (Kluwer Acad. Publ., 1990); Coniglio 
A., di Liberto F., Monroy G. and Perrugi F., Phys. Rev. 
B 44 (1991) 12605. 

A. Coniglio and H.J. Herrmann, Physica A, 225 1 (1996). 

M. Nicodemi and A. Coniglio, to be published. 

J. Arenzon, M. Nicodemi and M. Sellitto, to appear in J. 

Physique. 

For a recent review see Binder K. and Young P., Adv. in 
Physics 41 (1992) 547. 

J.B. Knight, C.G. Fandrich, C. Ning Lau, H.M. Jaeger, 
S.R. Nagel, Phys. Rev. E 51, 3957 (1995). 
A static configuration is defined by the criterion that 
during a fixed time t repo3e does not change any longer. 
We fixed in our simulation t repose — 330 time steps. Time 
is measured in such a way that unity corresponds to one 
single average update of all particles and all spins of the 
lattice. 

E. Ben-Nairn, J.B. Knight and E.R. Nowak, preprint in 



cond-mat/9603150 . 

M. Nicodemi, A. Coniglio, H.J. Herrmann, to be pub- 
lished. 

H. Hayakawa and D.C. Hong, proceedings of Nonlinear 
Dynamics and Chaos, Pohang, Korea, July 1995. 
T.A.J. Duke, G.C. Barker and A. Mehta, Europhys. Lett. 
13, 19 (1990). 

J. -P. Bouchaud, M.E. Cates and P. Claudin, J. Physique 
15, 639 (1995). 

C.-h. Liu, S.R. Nagel, D.A. Schecter, S.N. Coppersmith, 
S. Majumdar, O. Narayan, T.A. Witten, Science 269, 
513 (1995); S.N. Coppersmith, C.-h. Liu, S. Majumdar 



O. Narayan, T.A. Witten, preprint in :ond-mat/9511105. 
SF. Edwards and C.C. Mounfield, Physica A 226, 1 
(1996); ibid. 226, 12 (1996); ibid. 226, 25 (1996). 



3 



FIG. 1. Schematic picture of the lattice model considered 
here. Wavy and straight lines represent the two different kinds 
of bonds (eij = ±1). Filled (empty) circles are present parti- 
cles with spin Si = +1 (Si = —1). 

FIG. 2. Experimental data from [111 (square) and our 
MC data (circle) rescaled according eq. (H). Insert: density 
p(r, xo;t„) from our MC data as a function of tap number 
t n , for tap vibrations of amplitude xo = 0.001,0.01,0.05,0.1 
(from bottom to top) and duration r = 3.28 • 10 1 . The super- 
imposed curves are logarithmic fits from eq. (|). 

FIG. 3. The density profile p(h, r) as a function of depth 
h (h = corresponds to the top of the box, h — 200 to 
the bottom) for different values of the vibration duration r 
(r G [3.28-10- 3 ,4.92-10 4 ]). In the bulk of the system, for fixed 
h, p(h, t) is an increasing function of r. Continuous lines are 
Fermi-Dirac function fits from eq. (^). Insert: Rescaled den- 
sity profile p(h,r) / p s {r) as a function of the rescaled depth 
(h — ho(r))/s(T) (p s (t), ho(r) and s(t), arc fitting parameters 
to obtain the data collapse). 

FIG. 4. Force distribution P(v) as a function of weight v 
normalized by the mean force felt by the sites, for a static 
configuration of density p s — 0.764. Superimposed is the fit 
function in eq. (^). The fit parameters are a — 12 A, b — 5.6 
and c = 4.6. The distribution P(v) becomes narrower when 
the bulk density increases and is independent of the depth at 
which is measured (see [18]). 
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